﻿#include <algorithm>
#include <gtest/gtest.h>
#include <iostream>
#include <vector>

#include <xy/alg/integrator.h>

using namespace xy;
using namespace xy::math;
using namespace xy::alg;

TEST(TestIntegrator, NewtonCotes)
{
    {
        integrator_trapezoidal I(10000);

        auto f = [](double x) { return sin(numeric::pi * x); };
        EXPECT_FLOAT_EQ(I.solve(f, 0, 1), 2 / numeric::pi);
    }

    {
        integrator_simpson I(100);

        auto f = [](double x) { return sin(numeric::pi * x); };
        EXPECT_FLOAT_EQ(I.solve(f, 0, 1), 2 / numeric::pi);
    }

    {
        integrator_gauss I(5);
        // std::cout << I.m_nodes << std::endl;
        // std::cout << I.m_weights << std::endl;

        auto f = [](double x) { return sin(numeric::pi * x); };
        EXPECT_FLOAT_EQ(I.solve(f, 0, 1), 2 / numeric::pi);
    }

    {
        integrator_newton_tensor<3> I(5);
        auto f = [](const geom::vec<3, double> &x) { return math::vecxd{x[0] * x[1], x[1] * x[2], x[0] * x[2]}; };
        auto lower = geom::vec<3, double>{0.0, 1.0, 3.0};
        auto upper = geom::vec<3, double>{1.0, 3.0, 6.0};
        auto res = I.solve(f, lower, upper);
        EXPECT_LE((res - math::vecxd{6.0, 54.0, 13.5}).norm(), 1e-6);
    }
}